##

rm(list = ls())

library(tidyverse)
library(fixest)

## Helper functions 

source("code/helper_functions.R")

## Load Data

btw <- readRDS("data/data_main.RDS") %>%
  dplyr::select(
    county_id_2018,
    year,
    matches("zweit|diff"),
    east,
    treat_categorical4,
    matches("covar")
  ) %>%
  mutate(treat_categorical4 = factor(treat_categorical4,
    levels = c(
      "NoChange",
      "DecAndInc",
      "Increase",
      "Decrease"
    )
  ))

## Gen binary exit variable in t

btw <- btw %>%
  group_by(county_id_2018) %>%
  arrange(year) %>%
  filter(treat_categorical4 %in%
    c(
      "Decrease",
      "NoChange"
    )) %>%
  mutate(
    newspaper_exit =
      ifelse(treat_categorical4 == "Decrease", 1, 0)
  )



## Def covars

covars <- c("covar_pop_total", "covar_gdp_pc", "covar_emp_rate")
covars_diff <- paste0(covars, "_diff")

## Gen covar lags

btw <- btw %>%
  arrange(county_id_2018, year) %>%
  group_by(county_id_2018) %>%
  mutate_at(
    vars(all_of(covars)),
    list(lag = ~ lag(.))
  ) %>%
  ungroup()

## Covar lag list

cv_lags <- covars %>%
  paste0("_lag")

## Run models

cv <- cv_lags[1]

res <- lapply(cv_lags, function(cv) {

  ## Get DF

  df_temp <- btw

  ## Standardize covars

  df_temp <- df_temp %>%
    ungroup() %>% 
    mutate(across(one_of(cv), ~ . / sd(., na.rm = T)))

  ## Estimate

  f1 <- paste0(
    "newspaper_exit ~ ",
    cv,
    "| 0 "
  ) %>%
    as.formula()
  
  f2 <- paste0(
    "newspaper_exit ~ ",
    cv,
    "| year "
  ) %>%
    as.formula()
  
  f3 <- paste0(
    "newspaper_exit ~ ",
    cv,
    "| year +east "
  ) %>%
    as.formula()

  m1 <- feols(f1,
    data = df_temp,
    cluster = ~county_id_2018
  ) %>%
    tidy_feols() %>%
    mutate(fe = "none")
  
  m2 <- feols(f2,
    data = df_temp,
    cluster = ~county_id_2018
  ) %>%
    tidy_feols() %>%
    mutate(fe = "year")
  
  m3 <- feols(f3,
    data = df_temp,
    cluster = ~county_id_2018
  ) %>%
    tidy_feols() %>%
    mutate(fe = "year + east")

  ## Combine

  rbind(m1, m2, m3)
}) %>% reduce(rbind)

## Results

res_lev <- res %>%
  filter(!str_detect(term, "Intercept"))

## Same based on changes

nrow(btw)
btw %>%
  filter(!is.na(covar_emp_rate_diff)) %>%
  nrow()

cat("The mean of covar_emp_rate_diff is", mean(btw$covar_emp_rate_diff, na.rm = TRUE), "\n")
cat("The standard deviation of covar_emp_rate_diff is", sd(btw$covar_emp_rate_diff, na.rm = TRUE), "\n")
cat("The number of missing values in covar_emp_rate_diff is", sum(is.na(btw$covar_emp_rate_diff)), "\n")

cat("The mean of newspaper_exit is", mean(btw$newspaper_exit, na.rm = TRUE), "\n")

## Standardize covars

btw <- btw %>%
  ungroup() %>%
  mutate(across(
    all_of(covars_diff),
    ~ . / sd(., na.rm = T)
  ))

mean(btw$covar_emp_rate_diff, na.rm = T)

res_changes <- lapply(covars_diff, function(cv) {

  ## Get DF

  df_temp <- btw

  ## Estimate

  f1 <- paste0(
    "newspaper_exit ~ ",
    cv,
    "| 0"
  ) %>%
    as.formula()

  f2 <- paste0(
    "newspaper_exit ~ ",
    cv,
    "| year"
  ) %>%
    as.formula()

  f3 <- paste0(
    "newspaper_exit ~ ",
    cv,
    "| year + east"
  ) %>%
    as.formula()


  m1 <- feols(f1,
    data = df_temp,
    cluster = ~county_id_2018
  ) %>%
    tidy_feols() %>%
    mutate(fe = "none")

  m2 <- feols(f2,
    data = df_temp,
    cluster = ~county_id_2018
  ) %>%
    tidy_feols() %>%
    mutate(fe = "year")

  m3 <- feols(f3,
    data = df_temp,
    cluster = ~county_id_2018
  ) %>%
    tidy_feols() %>%
    mutate(fe = "year + east")

  ## Combine

  rbind(m1, m2, m3)
}) %>%
  reduce(rbind) %>%
  filter(!str_detect(term, "Intercept"))

res_changes %>%
  dplyr::select(estimate, term, fe) %>%
  filter(!str_detect(term, "diff_diff"))

## Remove all objects except for the results

rm(list = ls()[!ls() %in% c("res_lev", "res_changes")])

## Prep for plot

res_changes <- res_changes %>% 
  mutate(term = dplyr::recode(term,
    `covar_emp_rate_diff` = "Employment rate",
    `covar_gdp_pc_diff` = "GDP/capita",
    `covar_pop_total_diff` = "Population"
  )) %>%
  mutate(
    conf.low90 = estimate - qnorm(0.95) * std.error,
    conf.high90 = estimate + qnorm(0.95) * std.error
  ) %>%
  mutate(fe = dplyr::recode(fe,
    `none` = "Base model",
    `year` = "Year FE",
    `year + east` = "Year + East FE"
  )) %>%
  mutate(type = "Change between t-1 and t") %>%
    mutate(fe_order = case_when(
      fe == "Base model" ~ 1,
      fe == "Year FE" ~ 2,
      fe == "Year + East FE" ~ 3
    ))

## Results based on covar levels

res_lev <- res_lev %>%
  mutate(term = dplyr::recode(term,
    `covar_emp_rate_lag` = "Employment rate",
    `covar_gdp_pc_lag` = "GDP/capita",
    `covar_pop_total_lag` = "Population"
  )) %>%
  mutate(
    conf.low90 = estimate - qnorm(0.95) * std.error,
    conf.high90 = estimate + qnorm(0.95) * std.error
  ) %>%
  mutate(fe = dplyr::recode(fe,
    `none` = "Base model",
    `year` = "Year FE",
    `year + east` = "Year + East FE"
  )) %>%
  mutate(type = "Level in t-1") %>%
  mutate(fe_order = case_when(
    fe == "Base model" ~ 1,
    fe == "Year FE" ~ 2,
    fe == "Year + East FE" ~ 3
  ))

## Combine

res_full <- res_changes %>% 
  bind_rows(res_lev) %>% 
  mutate(type = fct_rev(type)) %>% 
  mutate(fe = fct_reorder(fe, fe_order))

## Plot this

pd <- position_dodge(0.4)

## Figure A16

p1 <- res_full %>% 
  ggplot(aes(term, estimate, type)) +
  geom_hline(yintercept = 0, linetype = 'dotted') +
  geom_errorbar(aes(ymin = conf.low90, ymax = conf.high90),
                width = 0, position = pd, linewidth = 1) +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high),
                width = 0, position = pd, linewidth = 0.5) +
  geom_point(shape = 21, position = pd,
             fill = 'white',
             size = 3) +
  facet_grid(fe~type, scales = 'free_x') +
  theme_bw()  +
  ylab('Effect on probability of newspaper exit')+
  xlab('Predictor')+
  coord_flip()
p1
